function [next_value] = contvalue( do_perturb, coeffs, fspace_next, index_node, t, max_time, statematrix, statematrix_next, cons_pc_node, weights_damage, weights_climsens, weights_cons, totalpop, statemin, Wiener, Params, option_ocean, option_damagetype)

% For "The Climate Risk Premium" (Lemoine, JAERE, 2020)

% Calculates continuation value at a given collocation node, for
% Epstein-Zin value function recursion

% coeffs contains coefficients for continuation value, to use with
% fspace_next

cons_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*cons_pc_node*Params.pop(t);
if Params.uselogstates==1
    M1_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*exp(statematrix(index_node,2));
    M2_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*exp(statematrix(index_node,3));
else
    M1_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,2);
    M2_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,3);
end
temp_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,4);

counter_fine = 1;

for index_interval = 1:1/Params.timestep
    
    % Calculate changes in variables
    index_time = t;
    sub_transitions
    
    % Update values of variables
    M1_value = M1_value + M1change;
    M2_value = M2_value + M2change;
    temp_value = temp_value + tempchange;
    cons_value = cons_value + conschange;
    
end % end stochastic process simulation

if do_perturb==1 % obtain consumption and temperature with exogenous increase in this period's CO2
    
    M1_store = M1_value;
    M2_store = M2_value;
    
    cons_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*cons_pc_node*Params.pop(t); % adjusts for state being per capita
    if Params.uselogstates==1
        M1_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*exp(statematrix(index_node,2)) + Params.perturbation*Params.M1_frac;
        M2_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*exp(statematrix(index_node,3)) + Params.perturbation*Params.M2_frac*exp(-Params.co2_decay*(t-1));
    else
        M1_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,2) + Params.perturbation*Params.M1_frac;
        M2_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,3) + Params.perturbation*Params.M2_frac*exp(-Params.co2_decay*(t-1));
    end
    temp_value = ones(size(Wiener.cons_increment,1),length(Params.warming),length(Params.cons_damage))*statematrix(index_node,4);
    
    counter_fine = 1;
    
    for index_interval = 1:1/Params.timestep
        
        % Calculate changes in variables
        index_time = t;
        sub_transitions
        
        % Update values of variables
        M1_value = M1_value + M1change;
        M2_value = M2_value + M2change;
        temp_value = temp_value + tempchange;
        cons_value = cons_value + conschange;
        
    end
    
    % restore CO2 trajectory
    M1_value = M1_store;
    M2_value = M2_store;
    
end

% store dimensions
sz = size(cons_value);
% obtain continuation value
if t < max_time-1
    % stack states into 1-D vectors
    cons_value = cons_value(:);
    M1_value = M1_value(:);
    M2_value = M2_value(:);
    temp_value = temp_value(:);
    % obtain matrix of states
    if Params.uselogstates==1
        next_state = [ log(cons_value/Params.pop(t+1)) log(M1_value) log(M2_value) temp_value ];
    else
        next_state = [ cons_value/Params.pop(t+1) M1_value M2_value temp_value ];
    end
    
    if Params.resetjumps==1
        temp_max = max(statematrix_next,[],1);
        temp_min = min(statematrix_next,[],1);
        for index_dim = 1:size(next_state,2)
            next_state(next_state(:,index_dim) > temp_max(index_dim),index_dim) = temp_max(index_dim);
            next_state(next_state(:,index_dim) < temp_min(index_dim),index_dim) = temp_min(index_dim);
        end
        %clear temp_max temp_min index_dim;
    end
    %     if sum( sum( next_state < min(statematrix_next,[],1)  ,1) ) > 0
    %         disp(['Jumped out below at state: ' mat2str(statematrix(index_node,:)) ', Min of next states: ' mat2str(min(next_state,[],1)) ', Min of next state grid: ' mat2str(min(statematrix_next,[],1)) ]);
    %     end
    % evaluate stored value function at those states
    next_value = funeval(coeffs',fspace_next,next_state);
    min_value = 1e-8*((1-exp(-Params.discount))*(Params.pop(t)/totalpop)*min(exp(statemin(:,1))).^(1-Params.time)).^(1/(1-Params.time));
    next_value = max( next_value, min_value ); % avoids negative values from extrapolating beyond last node
else
    cons_value = cons_value(:);
    next_value = ( (1-exp(-Params.discount))*(Params.pop(t+1)/totalpop)*(cons_value/Params.pop(t+1)).^(1-Params.time) ).^(1/(1-Params.time));    
end

next_value = next_value.^(1-Params.risk);
if Params.learn_ez~=1
    % return continuation value to previous dimensions
    next_value = reshape(next_value,sz);
    % take expectations over damage nodes
    next_value = bsxfun(@times,next_value,weights_damage);
    next_value = sum(next_value,3);
    % take expectations over climate sensitivity nodes
    next_value = bsxfun(@times,next_value,weights_climsens);
    next_value = sum(next_value,2);
end
% take expectations over consumption volatility
next_value = weights_cons'*next_value;

if ~isreal(next_value)
    error('next_value is complex')
end

end